functions {

	real GEVcpdf_lpdf(real x, real mu, real sig, real xsi){
		real z=(x-mu)/sig;
		real tau=exp(-6.0*abs(xsi));
		real lht;		
		if(1+xsi*z<tau){
			lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
		}
		else{
			lht=-log1p(xsi*z)/xsi;
		}
		return (1+xsi)*lht-log(sig);
	}
	
	real GEVtau_lpmf(int k, real x, real mu, real sig, real xsi){
		real z=(x-mu)/sig;
		real tau=exp(-6.0*abs(xsi));
		real lht;
		if(1+xsi*z<tau){
			lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
		}
		else{
			lht=-log1p(xsi*z)/xsi;
		}
		return -exp(lht)-lgamma(k+1);
	}
}

data {
	real xi_min;
    real xi_max;
	int<lower=0> T;
	array[T] int<lower=0, upper=6> nobs;
    vector[T] tau;
	matrix[6,T] y;   
    real sg_xi;
    real sg_alpha;
    real sg_s;    
}

parameters {
	// Level values for parameters
    real trans_xi_level;
    real ln_s_level;
    real m_level;
    
    // Innovation standard deviations (unscaled) for random walks
    real<lower=0.0,upper=5.0> g_xi;
    real<lower=0.0,upper=5.0> g_alpha;
	real<lower=0.0,upper=5.0> g_s;
    real<lower=0.0,upper=5.0> g_m;
    
    // Random walks
    vector [T] h_xi;
	vector [T] h_alpha;
    vector [T] h_s;
    vector [T] h_m;
}

transformed parameters {
    // set sg_m
    real sg_m = exp(ln_s_level);  // This preserves invariance to scale

    // Scale factors for innovations
    real gamma_xi = sg_xi*g_xi/sqrt(T);
    real gamma_alpha = sg_alpha*g_alpha/sqrt(T); 
    real gamma_s = sg_s*g_s/sqrt(T);
    real gamma_m = sg_m*g_m/sqrt(T);

     // vector of zeros
    vector [T] zeros = rep_vector(0.0,T);
    
    // Random walks
    vector [T] drw_xi = cumulative_sum(gamma_xi*h_xi)-mean(cumulative_sum(gamma_xi*h_xi));
    vector [T] drw_alpha = cumulative_sum(gamma_alpha*h_alpha)-mean(cumulative_sum(gamma_alpha*h_alpha));
    vector [T] drw_s = cumulative_sum(gamma_s*h_s)-mean(cumulative_sum(gamma_s*h_s));
    vector [T] drw_m = cumulative_sum(gamma_m*h_m)-mean(cumulative_sum(gamma_m*h_m));

   //  Construct Basic Parameters
    vector [T] trans_xi_1 = trans_xi_level + drw_xi;  
    vector [T] ln_alpha_1 = zeros;
    vector [T] ln_s_1 = ln_s_level + drw_s;
    vector [T] m_1 = m_level + drw_m;

    vector [T] trans_xi_2 = trans_xi_level + drw_xi;  
    vector [T] ln_alpha_2 = drw_alpha;  // note that level of ln_alpha is fixed at 0
    vector [T] ln_s_2 = ln_s_level + drw_s;
    vector [T] m_2 = m_level + drw_m;

    // Transformed Parameters
    vector [T] alpha_1 = exp(ln_alpha_1);
    vector [T] s_1 = exp(ln_s_1);

    vector [T] alpha_2 = exp(ln_alpha_2);
    vector [T] s_2 = exp(ln_s_2);
    
    // Construct GEV Parameters
    vector [T] xi_1 = xi_min + (xi_max-xi_min)*((exp(trans_xi_1))./(1+exp(trans_xi_1)));
    vector [T] sigma_1 = s_1.*(alpha_1.^xi_1);    
    vector [T] mu_1 = m_1 + (s_1./xi_1).*((alpha_1.^xi_1)-1);

    vector [T] xi_2 = xi_min + (xi_max-xi_min)*((exp(trans_xi_2))./(1+exp(trans_xi_2)));
    vector [T] sigma_2 = s_2.*(alpha_2.^xi_2);    
    vector [T] mu_2 = m_2 + (s_2./xi_2).*((alpha_2.^xi_2)-1);

}

model {
	g_xi ~ exponential(1.0);
    g_alpha ~ exponential(1.0);
    g_s ~ exponential(1.0);
    g_m ~ exponential(1.0);
    
    trans_xi_level ~ normal(0.0,1.5);
    // ln_s_level ~ normal(0.0,20.0);
    // m_level ~ normal(0,20.0);
	
    h_xi ~ std_normal();
	h_alpha ~ std_normal();
    h_s ~ std_normal();
    h_m ~ std_normal();
	
	for (t in 1:T) {
		// nobs[t] ~ GEVtau(tau[t],mu_1[t],sigma_1[t],xi_1[t]);
        nobs[t] ~ GEVtau(tau[t],mu_2[t],sigma_2[t],xi_2[t]);
		for (j in 1:nobs[t]) {
			// y[j,t] ~ GEVcpdf(mu_1[t],sigma_1[t],xi_1[t]);
            y[j,t] ~ GEVcpdf(mu_2[t],sigma_2[t],xi_2[t]);
		}
	}
}
